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Abstract. Multidimensional shock-capturing numerical schemes for special relativistic hydrodynamics (RHD) are 
computationally more expensive than their correspondent Euler versions, due to the nonlinear relations between 
conservative and primitive variables and to the consequent complexity of the Jacobian matrices (needed for the 
spectral decomposition in most of the approximate Riemann solvers of common use) . Here an efficient and easy-to- 
implement three-dimensional (3-D) shock-capturing scheme for ideal RHD is presented. Based on the algorithms 
developed by P. Londrillo and L. Del Zanna (Astrophys. J. 530, 508-524, 2000) for the non-relativistic magnetohy- 
drodynamic (MHD) case, and having in mind its relativistic MHD extension (to appear in a forthcoming paper), 
the scheme uses high order (third) Convex Essentially Non-Oscillatory (CENO) finite difference interpolation 
routines and central-type averaged Riemann solvers, which do not make use of time-consuming characteristic 
decomposition. The scheme is very efficient and robust, and it gives results comparable to those obtained with 
more sophisticated algorithms, even in ultrarelativistic multidimensional test problems. 
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1. Introduction 

Relativistic flows and shocks play an essential role in mod- 
ern high energy astrophysics, both for the interpretation 
of various observed features and for the description of the 
physical processes they give rise to, as, for example, the 
acceleration of highly energetic particles. Among the var- 
ious astrophysical objects in which relativistic flows have 
been invoked to explain the observed properties, the best 
studied are probably: 

1. Active galactic nuclei (AGNs), for which the presence 
of relativistic bulk motions (up to Lorentz factors of 
order 10) was soon suggested (Rees 1967). Associated 
to AGNs are often highly collimated relativistic jets 
(Begelman et al. 1984; and, for a review, Ferrari 1998), 
seen as apparent super-luminal lobes in some radio- 
loud sources and whose shock fronts are among the 
most promising candidates for high-energy particle ac- 
celeration. Mildly relativistic jets appear to be also as- 
sociated to galactic X-ray compact sources (generally 
called microquasars, see Mirabel & Rodriguez 1999, for 
a review). 
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Gamma-ray bursts (GRBs), supposed to be originated 
during the collapse of the iron core of massive stars 
and the subsequent fireball explosion (Piran 1999; 
Kobayashi et al. 1999), which give rise to an expanding 
blast wave of pairs and hadrons with typical Lorentz 
factors of 10 2 - 10 3 (Meszaros & Rees 1992), whose 
kinetic energy is then believed to be converted into 
gamma rays via cyclotron radiation and/or inverse 
Compton scattering. 

Pulsar wind nebulae, assumed to be bubbles of rel- 
ativistic particles and magnetic fields emitted by a 
pulsar as a relativistic wind with Lorentz factor rang- 
ing from 10 4 to 10 7 , depending on the model of pair 
production in the pulsar magnetosphere (Michel & Li 
1999). The wind region may be confined by a termi- 
nation shock, generated by the interaction with outer 
supernova matter, as is the case for the synchrotron 
emitting nebulae called plerions (Rees & Gunn 1974; 
Kennel & Coroniti 1984); or, if the pulsar is mov- 
ing with supersonic speed through the interstellar 
medium, by the resulting ram pressure, in this case 
giving rise to a bow-shock pulsar wind nebula that may 
be detected in Ha (Bucciantini & Bandiera 2001). 
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As we can see from these few but important examples, 
there is a strong interest in the astrophysical community 
to the development of computational codes for the numer- 
ical modeling and simulation of relativistic flows. 

Over the last decade, high resolution shock-capturing 
methods of Godunov type, successfully applied in classical 
fluid dynamics, have started to be employed for the case of 
relativistic hydrodynamics as well (Marquina et al. 1992; 
Schneider et al. 1993; Balsara 1994; Duncan & Hughes 
1994; Eulderink & Mellema 1994; Font ct al. 1994; Dolezal 
& Wong 1995; Falle & Komissarov 1996; Donat et al. 1998; 
Aloy et al. 1999). These schemes are characterized by the 
following main features: a conservative form of the dis- 
cretized equations, in order to capture weak solutions and 
satisfy jump relations; a reconstruction phase, to recover 
variables at inter-cell locations where fluxes have to be 
computed; an upwind phase, in which an exact or approx- 
imate solution to the local Ricmann problem is found. The 
simulated flows achieve high accuracy in smooth regions 
and, at the same time, shock profiles are stable and sharply 
defined. This is the reason for the great success of upwind 
differencing over central (or spectral) differencing. In the 
latter case, artificial viscosity must be introduced in order 
to damp the spurious oscillations that always form near 
discontinuities (Gibbs phenomena), thus leading to artifi- 
cial heating and unwanted damping of physical waves. 

When discontinuous solutions are of main interest, sec- 
ond order (both in time and space) total variation dimin- 
ishing (TVD) schemes coupled with an accurate Riemann 
solver are probably the best choice, in terms of sharp res- 
olution of discontinuities and overall stability. However, in 
the general case smooth wave-like or even turbulent fields 
appear together with discontinuous solutions, thus second 
order schemes are no longer able to resolve both features 
with enough precision. By relaxing the stringent TVD con- 
dition, a class of essentially non-oscillatory (ENO) higher- 
order schemes were first proposed by Harten et al. (1987), 
later modified in a more efficient implementation by Shu & 
Osher (1989). The ENO philosophy is to use adaptive sten- 
cils to reconstruct variables and fluxes at cell interfaces. 
Thus, in smooth regions symmetric stencils will be used, 
whereas near discontinuities the stencil will shift to the left 
or to the right, selecting the smoother part of the flow and 
thus achieving the same high resolution (typically from 
third to fifth order in the case of weighted ENO schemes, 
see Jiang & Shu 1996) everywhere, without the need of 
expensive adaptive mesh refinement (AMR) techniques. 
The price to pay is the presence, near shocks, of small 
oscillations of the order of the truncation error (Gibbs 
oscillations given by centered stencils that cross a discon- 
tinuity would be much higher, since they are proportional 
to the jump itself, independently on the resolution). 

Like for most TVD Godunov-type schemes, also for 
ENO algorithms the local spectral decomposition in the 
building of numerical fluxes is commonly adopted. It was 
strongly advocated already in the original paper, where 
a version of the ENO scheme with the (much) simpler 
component- wise reconstructions was found to fail, giv- 



ing significant oscillations, in some Riemann problems. 
However, due to the high resolution of ENO methods, 
characteristic decomposition on every point of the interpo- 
lation stencil becomes prohibitive when moving to three- 
dimensional simulations, especially for relativistic flows 
where Jacobian matrices are more complex than in the 
Eulerian case. The same kind of problem has to be faced 
in magnetohydrodynamics (MHD, see Londrillo & Del 
Zanna 2000, from now on LD), because of the increas- 
ing number of variables, equations and eigenmodes. The 
worst possible case is obviously that of relativistic MHD, 
for which only second order 1-D (Balsara 2001) and 2-D 
(Komissarov 1999) Godunov schemes have been presented 
so far. 

Following the scheme proposed in LD, in the present 
series of two papers a shock-capturing scheme that avoids 
the expensive characteristic decomposition, still retaining 
non-oscillatory properties, will be suggested, first for ideal 
relativistic hydrodynamics (RHD, in this paper) and then 
for ideal relativistic magnetohydrodynamics (RMHD, in 
the next paper of the series). Our method is based on the 
so-called central schemes, that extend the first-order local 
Lax-Friedrichs (LLF) flux splitting and assume an average 
over the Riemann fan at every cell interface (see Nessyahu 
& Tadmor 1990, for the original paper and Kurganov et 
al. 2001, for the latest developments in this field). The 
basic properties of our method can be found in Liu & 
Osher (1998), where a third order multidimensional cen- 
tral scheme was presented, with the following main fea- 
tures: 

1. Point values rather than cell averages are used (that 
is to say finite differences instead of finite volumes, see 
Sect. 3 for details), thus removing the need of staggered 
grids (employed in previous central schemes) and mak- 
ing the extension to multiple dimensions a trivial task. 

2. The semi-discrete form of the equations is solved, so 
that time integration can be achieved with any solver 
for ordinary differential equations (ODEs), for example 
TVD Runge-Kutta methods. 

3. No characteristic decomposition and Riemann solvers 
are required: fluxes are reconstructed and derived 
component- wise, thus achieving a great simplicity in 
the programming and, above all, efficiency. The only 
spectral pieces of information needed are the local 
highest characteristic speeds (related to the Courant 
coefficient). 

4. A new third order reconstruction algorithm is intro- 
duced and tested, called Convex-ENO (CENO), which 
has the fundamental property to stay as close as possi- 
ble to the lower order TVD limited linear reconstruc- 
tion (by using the minmod limiter in this case), thus 
reducing even to first order where needed, while re- 
taining high accuracy in smooth region. This was rec- 
ognized to be the key point of the surprising success of 
central schemes (Tadmor, cited in their acknowledge- 
ments), confirmed for the MHD case by LD. 
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Here, like in LD, we take advantage of all these positive 
features and we modify the original scheme by splitting the 
CENO reconstruction into two separate routines, so that 
reconstruction can be applied to primitive variables rather 
than to fluxes (the resulting scheme appears to be more 
robust, especially in the relativistic case). Moreover, dif- 
ferent limiters (e.g. the monotonized centered) and solvers 
(the upwind HLL, from Harten, Lax, and van Leer) are 
introduced and tested in our framework. 

We will show that our scheme, by providing an over- 
all high resolution, is able to compensate for the larger 
smearing of discontinuities, especially contact discontinu- 
ities (due to the use of solvers based on just one or two 
characteristic speeds) and the results are comparable with 
those obtained by more sophisticated (but much more 
computationally expensive, especially in 3-D) relativistic 
codes. Concluding, notice that the proposed scheme can 
be applied unchanged to Euler HD too, where only the 
definition of conservative variables, fluxes and character- 
istic speeds are different. 

2. Ideal RHD equations 

The covariant equations in special relativistic hydrody- 
namics (RHD) are (e.g. Landau & Lifshitz 1959): 



d a (pu a ) = 0, 
d a {wu"u p + pg aP ) = 0, 



(1) 

(2) 



where velocities are normalized against the speed of light 
(c — 1) and Greek letters indicate four- vectors, while Latin 
indexes will be devoted to spatial 3-D vectors. The metric 
tensor defines space-time properties and here a Minkowski 
flat space with g a/3 = diag{ — 1, 1, 1, 1} will be assumed 
throughout for simplicity, with coordinates x a = (t,x J ). 
The physical quantities involved in the conservation laws 
are the rest mass density p, the kinetic pressure p, the 
relativistic enthalpy w = e + p (e is the energy per unit 
volume, including rest mass energy) and the four-velocity 
u a = (7,71^'), where 7 = (1 — v 2 )^ 1 / 2 is the Lorentz 
factor. The system has to be closed with an equation of 
state p = p(p, e), and here the relation for an ideal gas 



p = (V - l)(e - p) => e = p + p/(T - 1) 



(3) 



will be considered, where the adiabatic index should be 
taken as L = 5/3 for the mildly relativistic case and as 
L = 4/3 for the ultrarelativistic case e>p. 

Godunov-type shock-capturing numerical methods de- 
veloped for classical Euler equations can actually be ap- 
plied to any multidimensional system of hyperbolic con- 
servation laws of the form 



du 
~dt 



d 

£ 

i=l 



dP(u) 

dx % 



(4) 



where u is the vector of conserved variables and /* are 
their corresponding fluxes, along each direction (d is the 



number of spatial dimensions). Equations (|l|) and 
automatically cast in this form by just defining 

u(v) = [p7 , W7 2 V , wj 2 — p] T , 



f(v) = [p"fv l , wj 2 vV + pS^ , u>7 



2„.ilT 



are 

(5) 
(6) 



where v = [p, , p) T are called primitive variables, and 
therefore the numerical techniques largely used for Euler 
HD equations can be applied to RHD too. The hyperbolic 
nature of Eq. (0) is guaranteed provided the local sound 
velocity c s is sub-luminal, that is for causal equations of 
state (Anile 1989). For F-law gases where Eq. (||) holds, 
the relativistic sound speed is given by 



dp 



a 2 I 1 



F - 1 



Tp/w, (7) 



which is obviously always less than unity, where s ~ pp~ r 
is the specific entropy and a 2 = Tpp -1 defines the classical 
sound speed. 

In any numerical time advancing routine, primitive 
variables v have to be derived from the conservative ones 
at least once per time step, and if this is trivial for Euler 
equations it is not so in the RHD case, in which a numeri- 
cal nonlinear root-finding technique must be employed. If 
we define W = wj 2 and 



u= [D,Q\E\ A 



(8) 



where D = p-f, Qi = Wv^ , and E = W — p, the system to 
be inverted can be cast in the single equation for 7: 



W(7) 2 (l-7~ 2 )-Q 2 = 0, 



(9) 



with w = p+T lP =>W= L>7+ri 7 2 p, here Ti = T/(r-l), 
and p ~ W — E, so that 



W{ 1 ) = 



£Fi7 2 - D7 
Ti7 2 -1 ■ 



(10) 



Once the Lorentz factor is found numerically with the re- 
quested accuracy, the primitive variables are easily recov- 
ered thanks to the relations above. 

3. A novel ENO-based central scheme 

In this section ENO methods and their implementation for 
component-wise central schemes, together with our spe- 
cific modifications, will be presented. For a general intro- 
duction and review of ENO methods for hyperbolic con- 
servation laws see Shu (1997). 

Consider a numerical discretization of Eq. (Q), in the 
one-dimensional case d = 1 to begin with. Given an inter- 
val [a, 6], N numerical cells U = [xi-1/2, 2^+1/2] of equal 
length Aa; = (b — a)/N can be defined, with cell centers 
(grid nodes) given by 



Xi = a + (i — 1/2) Ax; i = 1, . 



, N. 



(11) 



To an order r of spatial accuracy, the numerical value of 
any quantity v(x) (at a given time) will be denoted as 
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Vi = v(xi) + 0((Ax) r ) at grid points (the so-called point 
values) or Wj±i/2 = v(xi±i/ 2 ) + 0((Ax) r ) at cell bound- 
aries, where the order of accuracy refers only to smooth 
regions, in which the larger stencil can be used for inter- 
polation. Moreover, cell averages are defined as 



1 



Ax 



Vi = — / v(x)dx, 



(12) 



-1/2 



and only for schemes up to second order, r < 2, they do 
coincide with point values Vi. 

Most of the shock-capturing schemes evolve the cell 
averaged conserved quantities Ui in time, as obtained by 
integrating Eq. over the cell Ii. However, in the mul- 
tidimensional case, the resulting numerical fluxes / i; that 
discretize the physical ones along one direction, have to 
be averaged along the transverse directions. This implies 
a truly multidimensional numerical interpolation for high 
order schemes, which is usually complex and computa- 
tionally expensive. This approach is generally referred to 
as the finite volume approximation. Like in LD and in 
many other ENO schemes, from the works by Shu & Osher 
(1988; 1989) onwards, we will adopt here finite differ- 
ence approximations, based on point values. In the semi- 
discrete formalism (that is retaining continuous time de- 
pendency in the spatially discretized quantities), Eq. ([IJ) 
becomes (d = 1): 



du. t 
dt 



f 



i+l/2 



-1/2 



A;> 



(13) 



where f i±i/ 2 are high-order approximations to the primi- 
tives of physical fluxes, that is to say that cell averages of 
the f{x) function must coincide with point values f i of the 
flux function f{x), to the given accuracy (fi±i/ 2 = /i±i/2 
up to second order). The extension to the multidimen- 
sional case is now straightforward, since interpolations 
from cell centers to cell boundaries are made separately, 
dimension by dimension, and the other derivatives are just 
subtracted from the right hand side of Eq. (|l^) exactly in 
the same fashion as for the discretized x derivative. 

As specified by Shu & Osher (1988), Eq. @ must 
be integrated in time by using proper multi-level Runge- 
Kutta methods corresponding to the high order of spatial 
accuracy, so here we will always employ the optimal third 
order TVD algorithm: 



u 



n+1 _ 



u n + AtC[u n ] 

-u n + -u^ + -AtC[u^} 

4 4 4 

\u- + 2 -u^ + 2 -AtC[u^] 



(14) 



where the superscript n indicates the time stepping dis- 
cretization, it' 1 ) and refer to intermediate integration 
stages, and C[u] is simply the right hand side of Eq. (|l3|). 
The explicit time advancing scheme above is stable un- 
der the CFL (Courant-Friedrichs-Lewy) condition c < 1, 



where the Courant coefficient appears in the definition of 
the maximum time step allowed: 



At = 



(15) 



with a 1 being the largest speed of propagation of charac- 
teristic waves in the direction i. In the relativistic case, 
there is clearly a lower limit for the time step, given in 
Eq. ( |i~5| ) by simply taking a 1 = 1 along all directions. 

In order to complete the description of the scheme, the 
interpolation techniques and the approximate Riemann 
solver that defines the numerical inter-cell fluxes / i+1 / 2 
must be given. The following steps are taken, for every 
Runge-Kutta sub-cycle and for every direction in the mul- 
tidimensional case: 

1. Primitive variables are recovered from the conservative 
ones, according to the recipe already given in the pre- 
vious section (this step is actually taken just once for 
each sub-cycle): 



{ui} — ► {vi}; 



1, 



,N. 



(16) 



2. Primitive variables are reconstructed at cell interfaces 
to give two states, called left and right inter-cell states: 



{«<} 



{<i /2 }, {^f + i /2 }; 



o,. 



, N. 



(17) 



The two reconstructions are based on polynomial in- 
terpolation over different sets of stencils, one centered 
on Xi to approximate the left state vf +1 ^ 2 and the 
other centered on Xi+i to approximate the right state 
v H-i/2~ The interpolation is performed separately on 
each variable v(x) according to a Rec routine based 
on the CENO (Convex ENO, Liu & Osher 1998) tech- 
nique, as described in Appendix A. Contrary to com- 
mon ENO schemes, our method uses point values Vi to 
yield point value reconstructed quantities Wi+1/2, in- 
stead of starting from cell averages Cj. The order of 
the reconstruction is r = 3 in smooth regions, but it 
reduces to linear reconstruction or even to first order 
(by using minmod-type limiters) near discontinuities. 
In this latter case vf +1 f 2 = Vi an( ^ ^^+1/2 = f° r 
each discontinuous field. 
3. At each inter-cell point x i+ i/ 2 the local Riemann prob- 
lem must be solved, in some approximate way. The 
most accurate solution would be given by an exact 
solver that computes the evolved state from v L and 
v R and then the correspondent flux, to be finally iden- 
tified with /i+i/2- Another possibility is to define an 
average state v i+ i/ 2 and to decompose fluxes accord- 
ing to a linearized problem with the Jacobian df/du 
calculated there (Roe matrix approach). Our choice 
is to avoid spectral decomposition, which is compu- 
tationally expensive, and to take an average over the 
local Riemann fan as in central- type schemes, already 
mentioned in the introduction. A simple two-speeds 
Riemann solver is the HLL one, first used in central 
schemes by Kurganov et al. (2001), that retains an 
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upwind nature in the sense that coincides with f L or 
f R if the Riemann fan does not cross the inter-cell it- 
self (all the eigenvalues have the same sign). This may 
be written in the form 



pHLL 



a+f L + a-f R - a+a-{u R - u L ) 



(18) 



where all quantities are calculated from the recon- 
structed values of Eq. ( |l7| ) by using Eqs. (§) and (g). 
Here the or 1 coefficients take into account the highest 
speeds at the two sides of the Riemann fan, which can 
be estimated from the maximum and minimum eigen- 
value A of the Jacobians at the left and right states: 



a* = m&x{0,±\ ± (v L ),±\ ± (v R )}. 



(19) 



For relativistic flows, the required eigenvalues are (af- 
ter splitting the velocity along and perpendicular to 
the direction of spatial derivation, according to the rel- 
ativistic rule for addition of velocity vectors): 



A± = V- — !! , (20) 

1 — v z cj. 

that reduce simply to A* = (v±c 3 ) / (l±vc 3 ) in the one- 
dimensional case. Note that the maximum and mini- 
mum eigenvalues at the two reconstructed states are 
the only spectral pieces of information required. In 
this way shocks are handled correctly, whereas con- 
tact discontinuities and shear waves, corresponding to 
intermediate eigenvalues A = v\\ , can appear somehow 
smeared, when compared with the results from proper 
Riemann solvers. 

The simplest, smoothest, but also most dissipative nu- 
merical flux is the (local) Lax-Friedrichs one, given by 



cLLF 



f R - a(u R - u 1 



(21) 



in which a — max{a + ,a - } and therefore the averag- 
ing region is always symmetric with respect to 3^+1/2 
(the LLF numerical flux is the prototype of central 
schemes). Note that the maximum of the set of val- 
ues {cti+1/2} yields the CFL-related a coefficient ap- 
pearing in Eq. ([Hi]) , and actually they coincide for the 
so-called global LF numerical flux. The calculation of 
eigenvalues may be avoided completely in a global LF 
scheme that uses a = 1 everywhere, although this is 
the most smearing case. Note that Eq. ([l8]) reduces to 
Eq. ( plf ) when a + — aT = a, that is for symmetric 
Riemann fans in which vu = => A + = — A - . 
4. From the point values of numerical fluxes, the approx- 
imations of their derivatives must be finally calculated 
(for lower than third order schemes this step can be 
avoided) : 



{/i+1/2} > {/i+i/2} 



0,. 



,N. 



(22) 



These are the numerical fluxes that actually enter 
Eq. (p~3|) . This last step usually does not appear in 



other ENO high-order schemes, since the reconstruc- 
tion is made directly on fluxes and steps 2 and 4 
are taken simultaneously (the numerical flux is calcu- 
lated before reconstruction, using flux splitting meth- 
ods). However, we have noticed that for central- type 
schemes like the one presented here, which avoid char- 
acteristic decomposition, the proposed approach is 
more robust and less oscillatory. This final interpo- 
lation is again performed separately on each vari- 
able with another CENO-based routine (Der, see 
Appendix A). 

4. Numerical results 

In this section the proposed scheme is validated against 
typical tests available in the current literature, sepa- 
rated here in one-dimensional tests, essentially shock-tube 
Riemann problems, and a few 2-D and 3-D experiments. 
All the simulations have been run on a single PC-Linux 
(AMD Athlon 1GHz CPU, 512Mb RAM, Pacific-Sierra 
Vast/f90 compiler, -02 optimization and single precision 
calculations) in order to demonstrate that our scheme is 
not particularly demanding in terms of computational re- 
sources. In 3-D, with a grid of 100 3 nodes, typical speeds 
are of about 2 mins per Runge-Kutta iteration (three in- 
ternal sub-cycles for the third order scheme), so a sim- 
ulation like the spherical explosion takes a few hours of 
computational time. However, the code is fully parallelized 
with Message Passing Interface (MPI) directives and has 
been successfully tested on a variety of supercomputers, 
like CrayT3E, IBM SP3/4 and Beowulf clusters. 

In the following tests, we will always use the proposed 
third order CEN03 as a base scheme, unless specified oth- 
erwise. The approximate Riemann solvers employed are 
the HLL or the LLF, as described in the previous sec- 
tion, whereas the slope limiters used in the linear part of 
the reconstruction (see Appendix A) are the Monotonized 
Centered (MC) or the MinMod (MM). All the simulations 
presented here use the classical value T = 5/3 for the adi- 
abatic index. 

4.1. One-dimensional tests 

Shock-capturing numerical schemes, as in the classical hy- 
drodynamical case, must be able to reproduce the dis- 
continuous profiles involved in the so-called shock-tube 
problems. Given a unitary numerical 1-D pipe of N grid 
points, two constant states (p,v,p) are taken on the left 
(0 < x < 0.5) and on the right (0.5 < x < 1) with respect 
of a diaphragm, placed initially at x = 0.5 and then re- 
moved. Typical patterns seen in the subsequent evolution 
are shocks, contact discontinuities (characterized by den- 
sity jumps not accompanied by discontinuities in normal 
velocity and pressure) and rarefaction waves. In the rela- 
tivistic regime these features are qualitatively unchanged, 
since the structure of the characteristics is the same, but 
density jumps are not limited by any function of the adi- 
abatic index and rarefaction waves do not yield straight 
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Fig. 1. The relativistic blast wave problem 1 for time t = 
0.4. The computed profiles of density (diamond), pressure 
(cross), and velocity (plus) are shown against distance. 
The base third order CEN03 scheme is employed with 
MC limiter and HLL solver. The computational grid is 
formed by N = 400 equidistant cells and the Courant 
number is c = 0.5. 




Fig. 3. The same problem and settings as in Fig. 2, except 
for limiter and solver. Here the most smearing case of MM 
limiter and global LF flux splitting with a = 1 (highest 
possible value for relativistic flows) is tested. Differences 
arise only at x ~ 0.2 and x ~ 0.75, where the changes of 
slope due to the rarefaction wave are less sharp, and at 
the density peak, which is lower (p ~ 6.5) in this case. 



!last2: CEN03-HLL-MC, N = 400, CFL = 0.5 



1.2 


| i i i | i 

: p/1000 


i 


1 .0 


: i 










0.8 






%_ 


0.6 


\ / 
V / 


4> ~- 


0.4 




- 


0.2 


7 <°/ 10 / ^S^. 




0?f, 



left: 



Fig. 2. The relativistic blast wave problem 2 for time 
t = 0.35. The computed profiles of density (diamond), 
pressure (cross) , and velocity (plus) are shown against dis- 
tance. The base third order CEN03 scheme is employed 
with MC limiter and HLL solver. The computational grid 
is formed by N = 400 equidistant cells and the Courant 
number is c = 0.5. The exact value for the density peak 
is around 10.5, our numerical result is around 7.3, similar 
to what is found by other third order more sophisticated 
schemes. 



profiles, due to the nonlinear Lorentz transformation for- 
mulae. 

First we present two relativistic blast wave explosion 
problems, characterized by an initial static state with tem- 
perature and pressure much higher in the region on the 



Blast wave 1: 



Blast wave 2: 



(p,v, P ) L = (10,0,13.3), 
(p,v,p) R = (1,0, Hr 6 ), 

(p,v lP ) L = (1,0,1000), 
(p,v,p) R = (1,0,0.01), 



as in the Donat et al. (1998) paper. The first test is only 
mildly relativistic, while the second is more severe, with a 
shock speed corresponding to 7 ~ 6. In Fig. |l| and Fig. || 
the two tests are simulated with TV = 400 and Courant 
number c = 0.5, with the base scheme CEN03-HLL-MC. 
Note the total absence of oscillations and the accuracy in 
the definition of shocks and rarefaction waves. The contact 
discontinuity is more smeared, due to the use of a solver 
that does not take into account that intermediate wave 
and above all to the reconstruction method, which is the 
same for all the quantities, thus we cannot steepen (e.g. 
with the superbee limiter) the contact discontinuity alone, 
as usually done in characteristics based schemes. In par- 
ticular, the height of the density peak in Fig. || provides a 
measure of the numerical viscosity of the scheme: the ex- 
act value should be around 10.5, while here we find about 
7.3, slightly more than what shown in Donat's paper (in 
their third order simulation). A smaller value, and a larger 
spreading, is apparent in Fig. ^ where the most smear- 
ing case is tested, that is minmod limiter and global Lax- 
Friedrichs solver with a = 1, corresponding to A — > ±1 
(speed of light) in Eq. (|20|). However, given the extreme 
simplicity of the scheme, even this case should not be re- 
garded as completely unusable. 

The next example considered is the notorious relativis- 
tic shock reflection problem, where an ultrarelativistic cold 
wind hits a wall, a shock propagates backwards and a 
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Shock reflection: CEN02-HLL-MM, N = 250, CFL=0.5 Density perturbation: CEN03 — HLL — MC, N = 200, CFL = 0.5 
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Fig. 4. The relativistic shock reflection problem for time 
t = 0.75. The computed profiles of density (diamond), 
pressure (cross), and velocity (plus) are shown against 
distance. Here the second order CEN02 scheme is used 
with MM limiter and HLL solver. The computational grid 
is formed by N = 250 equidistant cells and the Courant 
number is c = 0.5. Note that the error in the density 
due to the wall heating phenomenon is around 2.3%, to 
be compared with the value of 2.5% given by Marquina's 
scheme in its third order implementation. 



static region of relativistically hot gas (e 3> p => <? s — > 
T — 1) is left behind. The numerical box is again [0, 1] (we 
use N = 250) and the reflecting wall is placed at x = 1. 
The physical values employed are: 

Shock reflection: (p, v,p) = (I, 0.99999, 0.01), 

that corresponds to a Lorentz factor as high as 7 ~ 223, 
about the highest allowed in single precision calculations. 
This is a very severe test, for the high velocities involved 
and for the so-called wall heating problem, visible as a 
dip in the density profile near the reflecting wall. This 
is a classical numerical artifact, due basically to the im- 
plicit numerical viscosity present in every scheme. In Fig. ^| 
a simulation with second order accuracy (both in time 
and space), HLL solver and MM limiter is shown for time 
t = 0.75 (c = 0.5). This is the only case where the third 
order CEN03 has failed, giving significant postshock oscil- 
lations, that can only be reduced by lowering the Courant 
number or by enhancing numerical viscosity, though they 
never completely disappear. However, the overheating er- 
ror in the density is here only 2.3%, in spite of the sec- 
ond order of accuracy and of our simplified scheme, to be 
compared with the 2.5% value of Donat's paper in their 
third order implementation of the celebrated Marquina's 
scheme. 

As a last 1-D application we present a test proposed 
in the Dolezal & Wong (1995) paper, which nicely shows 
the capacity of ENO-based schemes of treating accurately 
both discontinuous and smooth features occurring close 
together and at the same time. A shock tube is perturbed 




Fig. 5. The Riemann problem with sinusoidal density per- 
turbation for time t = 0.35. The diamonds are the results 
obtained with CEN03, MC limiter and HLL solver, in a 
simulation with N = 200 grid points and Courant number 
is c = 0.5. The solid line is the density profile that comes 
out from a simulation with 2000 grid points. High-order 
ENO schemes are particularly suited for problems involv- 
ing shocks interacting with smooth wave-like structures. 

in its right-hand state (0.5 < x < 1) with a density sinu- 
soidal profile: 



Density perturbation: 



(p,v,p) L = (5,0,50), 

(p,v,p) R =(2 + 0.3sin50x,0,5). 



The subsequent evolution in time shows the interaction 
of the blast wave with the density wave (the values are 
not exactly the same as in Dolezal's paper, where a nu- 
clear equation of state is used), which enters the expand- 
ing heated region and modulates its density plateau. In 
Fig. H the extreme accuracy of our third order CEN03 
scheme is shown, by comparing the resulting density pro- 
files at time t = 0.4 in two runs, one with just N — 200 
points (diamonds), tested against one with N = 2000 
grid points (solid line). Again, by comparing our results 
to those obtained by the third order characteristics based 
ENO scheme in Dolezal's paper, it is apparent that no sig- 
nificant differences arise, in spite of the much less effort 
involved in our scheme. 

4.2. Multidimensional tests 

Multidimensional relativistic simulations are more diffi- 
cult than one-dimensional ones because the velocity com- 
ponents are spatially interpolated separately, possibly 
causing the condition v 2 < 1 to fail in ultrarelativistic 
regimes due to numerical errors in the reconstruction. For 
this reason in some cases we had to reduce to first order 
reconstruction, namely when the Lorentz factor reaches 
7 = 10. Note that this does not mean in any sense that 
this threshold cannot be exceeded, but only that the accu- 
racy is lower in these regions, which usually are located at 
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Ricmann2D: CENOJ-LLF 



N = 400, CFL = 0.b 




Fig. 6. The relativistic 2-D Riemann problem for time 
t = 0.4. Density contours (30) in logarithmic scale are 
shown, for a simulation using CEN03, MM limiter and 



LLF solver, with N x 
number is c = 0.5. 



400 grid points and Courant 



shock fronts where the order would have been lowered any- 
way. We have also tested reconstruction on four-velocity 
components, which are not bounded, but in this way the 
problem is just shifted to the primitive variables finding 
routine, thus with no improvement at all. 

While it easy to test one-dimensional codes, since 
Riemann problems can be solved exactly through iterative 
algorithms, it is not so in more than one dimension, where 
it is rare to see really quantitative numerical scheme vali- 
dations. Here we will propose a 2-D Riemann problem, 2-D 
and 3-D explosions, compared with the correspondent 1-D 
cylindrical and spherical solutions, and finally one simu- 
lation of a relativistic jet, which is not a quantitative test 
but it is a true astrophysical application. 

The two-dimensional counterpart of shock tubes is a 
square domain divided in four quadrants of constant val- 
ues at initial time and free evolution for t > 0. In Lax & 
Liu (1998) all the possible different configurations involv- 
ing 1-D shocks, 1-D rarefactions waves and 2-D slip lines 
(contact discontinuities) were studied in detail. In Fig. ^| 
we show the output (contours of the density logarithm for 
time t = 0.4) of a situation similar to their configuration 
12, where the four boundaries defines two contact discon- 
tinuities and two 1-D shocks (symmetric with respect to 
the main diagonal) . Here a relativistic version of this test 



is proposed, with the following initial settings: 

(p,v x ,v y ,p) NE = (0.1,0,0,0.01), 

(p,v x ,v y ,p) NW = (0.1,0.99,0,1), 

(p,v x ,v y ,p) 3W = (0.5,0,0,1), 

(p,v x ,v y ,p) SE = (0.1,0,0.99,1). 



Riemann 2-D: 



Note that we have not taken exact 1-D shocks across the 
N and E interfaces, and this may be recognized by ob- 
serving the evolved discontinuities in Fig. q, converging 
towards the NE corner with their complete Riemann fans. 
In the rest of the domain the structure evolves with curved 
shock fronts and a complicated pattern in the SW quad- 
rant, reminiscent of an oblique jet with bow shock and 
converging internal shock fronts. The lines in the SW di- 
rection with respect to the bow shock are actually due to 
spurious waves created at the W and S interfaces by the 
numerical diffusion term in the energy equation (left and 
right states have a jump in kinetic energy), and cannot be 
removed if not by using a Roe- type solver. Note that the 
most diffusive case (LLF solver and MM limiter) was used 
in the simulation, since other cases resulted more unstable 
at the same high resolution and with the same Courant 
number c = 0.5. 

Fig. ^ and Fig. prefer respectively to the 2-D cylindri- 
cally symmetric and 3-D spherically symmetric blast wave 
explosion problem below: 



Radial blast wave: 



{p,v r ,p) = (1,0,1000) r<0.4, 
{p,v r ,p)= (1,0,1) r>0.4, 



where in both cases the runs have been performed in a 
Cartesian unit box (with reflective right hand side bound- 
ary conditions along each direction) and compared with 
the correspondent 1-D radial solution, obtained with N = 
800 grid points in cylindrical or spherical coordinates. The 
results are satisfactory, since the radial symmetry is well 
preserved (here the results along the diagonal are shown) 
and high Lorentz factors are reached without particular 
problems (7 > 25 in the spherically symmetric case). Here 
the scheme used is CEN03 with HLL solver and MC lim- 
iter, with c = 0.5 and c = 0.3 in the 2-D and 3-D cases, 
respectively. 

Finally, as a typical astrophysical test, we simulate the 
propagation of an axisymmetric jet in 2-D cylindrical co- 
ordinates (z, r). Note that jet simulations are a very hard 
test for codes not based on characteristics decomposition, 
because of usually stronger numerical viscosity at shear 
layers. The domain is < r < 8 and < z < 20, with re- 
flective boundary conditions on the axis r — and simple 
extrapolation at the other boundaries (except at z = 
within the jet radius, where initial values are kept con- 
stant). At t = a relativistic jet with v z — 0.99 and 
density 100 times less than the surroundings (but same 
pressure) is located at r < 1 and z < 1 



Jet: 



{p,v z ,v r ,p) = (0.1,0.99,0,0.01) r < l,z < 1, 
(p,v z ,v r ,p) = (10,0,0,0.01) outside; 



density and velocity jumps are actually smoothed in or- 
der to reduce spurious transverse waves that appear due 
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Blast2D: CEN03-HLL-MC, N = 250, CFL = 0.5 



Blast3D: CEN03-HLL-MC, N = 100, CFL=0.3 
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Fig. 7. The relativistic 2-D blast wave problem for time 
t = 0.4. Here the simulation employed CEN03 reconstruc- 
tion, MC limiter and HLL solver, with N x = N y = N = 
250 grid points. The domain is a Cartesian 2-D unit box 
[0,1] x [0,1], with initial enhanced pressure (a factor of 
10 3 ) in a disk sector centered at the origin with radius 
r max = 0.4 (reflecting boundary conditions are assumed 
for x — and y = 0). A radial cut (along the main diago- 
nal) of the computed quantities are compared with those 
obtained through a high-resolution (N = 800) 1-D sim- 
ulation in cylindrical coordinates (solid line), using the 
same parameters. For a more quantitative comparison, in 
the bottom panel the density relative error is plotted. Its 
rather large value near shocks is simply due to the natu- 
rally reduced accuracy at discontinuities and to the very 
small number of grid points (eight) present in the density 
peak. 



to numerical viscosity at the shear layer. The material is 
injected with a Lorentz factor 7 ~ 7.1, corresponding to 
a relativistic Mach number (e.g. Duncan & Hughes 1994) 
of M. — "fv/j Cv c s ~ 17.9. The jet evolution is followed 
until t = 40, as shown in Fig. ^, where density contours 
and gray shades in logarithmic scale are presented. The 
code settings are CEN03-HLL-MC with Courant coeffi- 
cient c = 0.25, while the resolution employed is 400 x 160, 
corresponding to 20 grid points per jet radius. Note that 
the smearing of contact discontinuities, unavoidable in 
methods not based on characteristics decomposition, is ac- 
tually small and vortices due to Kelvin-Hclmoltz instabil- 
ities are nicely defined, as well as the external bow shock, 
the internal Mach disk and other shocks reflected off the 
axis. Moreover, notice the absence of the so-called carbun- 



Fig. 8. The relativistic 3-D blast wave problem for time 
t = 0.4. All the settings are the same as in Fig. |?], except 
that here N = 100 along all directions and the Courant 
number is lower (c = 0.3). A radial cut (along the main 
diagonal) of the computed quantities are compared with 
those obtained through a high-resolution (N = 800) 1-D 
simulation in spherical coordinates (solid line), using the 
same parameters. 



cle problem, usually manifesting as an extended nose in 
the front of the jet on the axis (e.g. Quirk 1994). 

5. Conclusions 

The central shock-capturing scheme of Londrillo & Del 
Zanna (2000) is extended to the relativistic case, here 
for ideal RHD flows and to RMHD in the next paper of 
the series. Compared to other schemes proposed for rela- 
tivistic astrophysical problems over the last decade in the 
literature, our method is extremely simple and efficient, 
since no eigenvector decomposition and Riemann solvers 
are involved. In spite of this, due to the high order ac- 
curacy achieved in smooth regions and to TVD limiting 
near discontinuities, provided by the CENO reconstruc- 
tions employed, our results compete with those obtained 
by more sophisticated (but less computationally efficient) 
codes, even in 2-D and 3-D highly relativistic test prob- 
lems. Different geometries (cylindrical and spherical coor- 
dinates) and boundary conditions have also been tested. 
All the numerical experiments have been run on a simple 
PC-Linux machine, even the most demanding 3-D test 
with one million grid points, in order to demonstrate the 
efficiency of our code. We believe that our simple scheme 
can be successfully employed in many relativistic simula- 
tions of astrophysical interest, either in the fluid case and 
above all in the most demanding magnetic case, which will 
be the subject of a forthcoming paper. 
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Jet: CEN03-HLL-MC, N,=400, N = 160, CFL = 0.25, 




5 10 15 20 



Fig. 9. The relativistic jet for times t = 20, 30, 40 in 2- 
D cylindrical geometry (logarithmic density contours and 
shades are shown). The inflow speed is v z = 0.99 and the 
internal to external density ratio is r\ — f/fOO, with an 
overall pressure equilibrium, leading to a beam relativistic 
Mach number of 17.9. The jet radius is assumed as length 
unit, corresponding to 20 computational cells (N z = 400 
grid points in the axial direction and N r = 160 points in 
the radial direction). CEN03 reconstruction, HLL solver 
and MC limiter have been used in the simulation, with 
Courant number c = 0.25. 



Appendix A: Convex ENO interpolation routines 

Given a function v(x) with grid point values {u;}, for 
i = 1,...,N, let us see how non-oscillatory reconstruc- 
tions (Rec) and second derivatives (Der) are defined in 
our scheme. These procedures are based on the original 
CENO reconstruction from Liu & Osher (1998) and they 
have been described also in LD (we repeat them here for 
completeness) . 

In the range £j_i/2 < Xi < x i+1 / 2 , we first consider 
the two linear polynomials Li(x), based on the stencils 
[£j_i,a;j] or [a^, and we choose their convex com- 



bination which is closest to Vi, that happens to coincide 
with the TVD limited reconstruction: 

L i (x)=v i + v' i (^p-y (A.l) 

Here the non-oscillatory undivided first derivative v[ is 
given by the minmod function: 

vl =mm[A_» i) A + w j ], (A.2) 

and A±Vi = ±(«i±i ± 1). Thus, to second order of ac- 
curacy the CENO technique simply reduces to the usual 
TVD monotone slope limiting. Although not in the philos- 
ophy of the original CENO approach, we have tested other 
limiters. The most compressive ones, like superbee are too 
oscillatory in our component- wise framework, while a good 
compromise between sharpening and stability appears to 
be provided by the so-called monotonized centered (TVD 
for 1 < 9 < 2): 

vl = mmpA-Wi.AoVi.flA+Vj], (A.3) 

where A Vi — (A + Vi + A_Vi)/2 = (v i+ i — Vi-\)/2 and 
9 = 2 in order to bias towards central interpolation (since 
twice the one-sided derivatives might be chosen near dis- 
continuities, the Courant coefficient should be lowered cor- 
respondingly for overall stability). In both relations the 
function mm is defined as 

{minj{aj} aj > Vj, 
m&Xj{a,j} aj < Vj, (A. 4) 

otherwise. 

To achieve an accuracy of third order in smooth re- 
gions, we need to calculate the three quadratic polynomi- 
als (k = -1,0, +1, j = i + k): 

QH X )=v 3+ A aVj +I A fo (^)* (A.5) 

where A§«j = A + A_«j = A_A + Uj = Vj + i — 2vj + Vj-i, 
together with the corresponding weighted differences 

d\ = a k {Q k i -U). (A.6) 

Then, we take again their convex combination which is 
closest to Li at a given point x (this procedure is the 
generalization of minmod limiting to higher orders) and 
finally, in smooth regions where all d\ have the same sign, 
we take 

Qi = Q\\ |rfM=min fe (|^|), (A.7) 

otherwise we exit from the interpolation routine with 
Qi = Li. As for the coefficient 9 in the limiter of the lower 
order, the weights a k are chosen in order to bias towards 
central interpolation, and following the original recipe for 
component-wise CENO schemes we take a* 1 = 1 and 
a = 0.7. In our Rec routine, the left and right states are 
calculated at the same time and the output is finally 

v i+l/2 = Qi{Xi+\/2), V?-l/2 = Qi{Xi-l/2)- (A.8) 
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{vi} — ► Vi=Vi- ^v'/. (A.9) 

In this case, the above CENO selection criterion applies 
with null linear polynomials, = 0, and quadratic ones 
given by (k = — 1, 0, +1, j = i + k): 

Qi = A„ vj = v j+1 - 2vj + Vj-i, (A.10) 

then finally v" = Qi (which is zero, so Vi = Vi as for second 
order approximation, in the case different signs of the Q\ 
terms). 
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The other CENO routine used in our scheme is 
Der, which allows one to calculate non-oscillatory second 
derivatives v", needed to transform cell averages into point 
values, calculated at the same point: 



